#include "kernel.h"
#include "channeler.h"
#include "variables.h"
#include "data.h"
#include "psrfold.h"

void prepare(long block_size,int ntap,int nfft)
{
    pol1   = malloc(sizeof(double) * block_size);
    pol2   = malloc(sizeof(double) * block_size);
    opol1  = malloc(sizeof(double) * block_size * 2);
    opol2  = malloc(sizeof(double) * block_size * 2);
    coeffs = malloc(sizeof(double) * nfft * ntap);
    ocoeffs = malloc(sizeof(double) * nfft * (2 * ntap - 1));
    half_data1 = calloc(sizeof(double) , nfft / 2);
    half_data2 = calloc(sizeof(double) , nfft / 2);

    psize = get_period_size(400.0 / (1.0 * nfft / 2));
    // psize = get_period_size(400.0);

    profile_out = calloc(sizeof(double) , (nfft / 2 + 1) * nbin);

    subband1 = malloc(sizeof(fftw_complex) * block_size * 2);
    subband2 = malloc(sizeof(fftw_complex) * block_size * 2);

    pdata = calloc(sizeof(double),psize * (nfft / 2 + 1));
    pnum  = calloc(sizeof(double),psize * (nfft / 2 + 1));

    // pdata = calloc(sizeof(double),psize);
    // pnum  = calloc(sizeof(double),psize);

    location = calloc(sizeof(long),nfft / 2 + 1);
}


void operation(const char *input_name,const char *output_name,const char *coeffs_name,
               int ntap,int nfft,long block_size,int opfb_factor)
{
    prepare(block_size,ntap,nfft);
    FILE *input_pf = load_data(input_name);
    load_coeffs(coeffs_name,ntap,nfft,coeffs);
    if(opfb_factor == 2){
        /* oversampling pfb */
        opfb_coeffs(coeffs,ocoeffs,ntap,nfft);
        ntap = 2 * ntap - 1;
    }
    if(ntap * nfft> block_size){
        fprintf(stderr,"operation():the block size(%ld) is samller than ntap * nfft(%d)\n",block_size,ntap*nfft);
        exit(0);
    }
    printf("after oversample, nfft = %d,ntap = %d\n",nfft,ntap);
    unsigned nblock = (get_file_length(input_name) - 4096) / (2 * block_size);
    printf("nblock: %d\n",nblock);

    long res_num = 0;
    long bblock_size = block_size;
    nblock = 1;
    for(int i = 0;i < nblock;++i){
        printf("process %d blocks\n",i);
        unpack_dada(input_pf,pol1,pol2,block_size);
        if(opfb_factor == 2){
            opfb_data(pol1,opol1,half_data1,nfft,block_size);
            opfb_data(pol2,opol2,half_data2,nfft,block_size);
            bblock_size = block_size * 2;
            res_num = bblock_size / nfft;
        }
        if(opfb_factor == 1){
            bblock_size = block_size;
            res_num = bblock_size / nfft;
        }
        printf("begin channeler\n");
        channeler(opol1,ocoeffs,subband1,res_num,ntap,nfft);
        channeler(opol2,ocoeffs,subband2,res_num,ntap,nfft);   
        printf("begin dedispersion\n");
        dedispersion_fold(subband1,subband2,pdata,pnum,location,nfft,psize,bblock_size);
    }
    printf("begin integral_data\n");
    for(int i = 0;i < (nfft / 2 + 1);++i){
        if(i == 0 || i == (nfft / 2)){
            integral_data(pdata+i*psize, profile_out+i*nbin, psize / 2);
        }else{
            integral_data(pdata+i*psize, profile_out+i*nbin, psize);
        }
    }
    write_data(output_name,profile_out,(nfft / 2 + 1) * nbin);
    destroy();
}

void destroy()
{
    free(pol1);
    free(pol2);
    free(opol1);
    free(opol2);
    free(coeffs);
    free(ocoeffs);
    free(pdata);
    free(pnum);
    fftw_free(subband1);
    fftw_free(subband2);
    free(profile_out);
    free(half_data1);
    free(half_data2);
    free(location);
}